function vnsAnalysis(inDir,outDir)
%inDir='/Users/matthew/workspace/mfiles/projects/data/VNSrawdata/8sec/';
%outDir='/Users/matthew/workspace/mfiles/projects/data/VNSpli/8sec/';


pliXBand(inDir,outDir);


function [data]=import2FT(filename)
fsample=256;
trialLen=2048;
nTrials=3;
nch=19;
filename
toImport=importdata(filename,'\t',2); 

toImport.data=toImport.data';
data.label=toImport.colheaders(1:19);
data.fsample=fsample;
sample_interval=1/fsample;

%[nSamples,nch]=size(toImport.data);
%nTrials=nSamples/trialLen;
data.trial=cell(1,nTrials);

indx=1;
indx=1:trialLen:trialLen*nTrials;
for i=1:nTrials 
    data.trial{i}=toImport.data(1:nch,indx(i):indx(i)+trialLen-1);
    data.time{i}=(0:1:length(data.trial{i})-1)*sample_interval;
end


function [bf_trials]=selectBand(trials,band)

bf_trials=trials;
nTrials=length(trials.trial);
cfg = [];
if(band(1)~=0 && band(2)~=0)
    cfg.bpfilter='yes';
    cfg.bpfreq=band;
end    
bf_trials= ft_preprocessing(cfg,trials);

function [pli3d,spli3d]=pliXTrial(bf_trials)

nTrials=length(bf_trials.trial);
nChs=size(bf_trials.trial{1},1);
pli3d=zeros(nTrials,nChs,nChs);
spli3d=zeros(nTrials,nChs,nChs);
for i=1:nTrials
    [pli3d(i,:,:),zscorematrix,spli3d(i,:,:)]=pli2(bf_trials.trial{i}');  
end

function pliXBand(inDir,outDir)

bands={[1 4],[4 8],[8 13],[13 30],[30 49]};
bandLabels={'delta','theta','alpha','beta','gamma'};
nBands=length(bandLabels);
filter=fullfile(inDir,'*.dat');
files=dir(filter);
nSubjs=length(files);


for i=1:nSubjs
    subjName=files(i).name;
    trials=import2FT(fullfile(inDir,subjName))
    
    for j=1:nBands   
        
        bf_trials=selectBand(trials,bands{j});
        
        [subjName,dot]=strtok(subjName,'.');
     
        
        [pliXT ,spliXT]=pliXTrial(bf_trials);
        
        outfile=strcat('pli_',bandLabels{j});
        outfile=strcat(outfile,'_');
        
        outfile=strcat(outfile,subjName);
        
        save(fullfile(outDir,outfile),'pliXT','spliXT');
    end    
end


function [PLI,zscoreMatrix,sPLI]=pli2(m)

% a is a filtered multichannel signal (time x channels)
% hilbert(a) calculates analytic signal (complex valued) of each
% column of a. Phase Lag Index between channel i and j averaged over
% time bins is stored in PLI(i,j)
% number of channels

nChs=size(m,2);
timeStop=size(m,1);

PLI=zeros(nChs,nChs);
complex_m=hilbert(m);

phi=atan2(imag(complex_m),real(complex_m));


zscoreMatrix=zeros(nChs,nChs);

for i=1:nChs
    
    for j=1:nChs
        count=0;
        count2=0;

        for t=1:timeStop
        
            phaseDiff=phi(t,i)-phi(t,j);
            
            if(sin(phaseDiff)>0)
                count=count+1;
            end
            if(abs(sin(phaseDiff))>0)
                count2=count2+1;
            end
            
        end
        
        if(count2>0)
            pli=2*abs(.5-1*count/count2);
            %for FDR
            zscoreMatrix(i,j)=((count/count2)-0.5)/sqrt((0.5*0.5)/count2);
        else
            pli=0;
        end
        
        PLI(i,j)=pli;
    end
end
sPLI=fdr(zscoreMatrix,PLI,0.01);
function [pliM]=fdr(zscoreMatrix,pliM,alpha)

lindx=find(triu(ones(size(pliM)),1));
[row,col]=ind2sub(size(pliM),lindx);
z=abs(zscoreMatrix(lindx));
k=ones(1,length((z)));%zeros or ones?? stam one
N=length(z);

for i=1:length(z)
    
    k(i)=k(i)+sum(z(i)<z);%to be equal of stam code
    
    if(z2p(z(i))>(k(i)*alpha)/N)
        pliM(row(i),col(i))=0;
        pliM(col(i),row(i))=0;
    end        
        
end
function [p]=z2p(zscore)
        p=1;
       
        if(zscore > 0.67 )
            p=0.5;
        end
        if(zscore > 0.84)
            p=0.4;
        end
        if(zscore > 1.04 )
            p=0.3;
        end
        if(zscore > 1.28 )
            p=0.2;
        end
        if(zscore > 1.64)
            p=0.1;
        end
        if(zscore > 1.96 )
            p=0.05;
        end
        if(zscore > 2.33 )
            p=0.02;
        end
        if(zscore > 2.58 )
            p=0.01;
        end
        if(zscore > 2.81 )
            p= 0.005;
        end
        if(zscore > 3.09 )
            p=0.002;
        end
        if(zscore > 3.29 )
            p=0.001;
        end
        if(zscore > 3.89 )
            p=0.0001 ;
        end



